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ABSTRACT 

MHD simulations of sunspots have successfully reproduced many aspects of sunspot fine structure 
as consequence of magneto convection in inclined magnetic field. We study how global sunspot 
properties and penumbral fine structure depend on the magnetic top boundary condition as well 
as on grid spacing. The overall radial extent of the penumbra is subject to the magnetic top boundary 
condition. All other aspects of sunspot structure and penumbral fine structure are resolved at an 
acceptable level starting from a grid resolution of 48 [24] km (horizontal [vertical]). We find that 
the amount of inverse polarity flux and the overall amount of overturning convective motions in the 
penumbra are robust with regard to both, resolution and boundary conditions. At photospheric levels 
Evershed flow channels are strongly magnetized. We discuss in detail the relation between velocity 
and magnetic field structure in the photosphere and point out observational consequences. 
Subject headings: MHD - convection - radiative transfer - sunspots 



1. INTRODUCTION 

Numerical simulations of sunspot fine structure have 
seen substantial improvem ents over that past 5 years. 



Schiissler fc Vogler] ( 2006 ) presented the first realistic 
MHD simulations or umbral do ts, followed by a series o f 



simulations in ' slab-geo metry': Hcincmann ct al. (2007); 
Scharmer et al.| (]2008[); Rempel et al. (2009b) addrcss- 



ing the transition from umbra to penumbra and gave the 
first insight into the or igin of penumbra fine s tructure 



and the Evershed flow. Kitiashvili et al. (2009) focused 



on a magneto convection simulation in strongly inclined 
field representative of penumbral conditions and studied 
the influence of field stren gth and inclination on the pres 
ence horizo ntal outflows. 
( 2011a|b|c ) focused on ful 



Rempel et al. (2009a); Rempel 
circular sunspots with mcreas- 
ing realism, wi th the best reso lved simulations such as 
those shown in Rempel (2011a) reaching and exceeding 
the richness of detail seen in the best available obserya- 
tions today. As described in detail in Rempel ( |2011a|b l 
these simulations reproduce the key aspects 01 penum- 
bral fine structur e as inferred from high resolution obser- 
vations (see, e.g. Scharmer et al. 2002 Langhans et al.| 
20051 |RmmiekTfc Maxino|liMOfl| |lchimoto et ai.J[g(K)7a|b| 
Langhans et al.||2007| |Scharmer et al.||2 007 ; Rimmclc 
2008; Franz & Schlichenmaier 2009; Bcllot Rubio ct al. 
2010; Franz 2011) and recent reviews by Solanki (2003); 



Thomas & Weiss (2004, 2008); Rempel & Schlichenmaier 
( |201l[ ); jBorreroinclirmoto| |2011| ) 



Simulations are still tar from perfect. The numerical 
grid resolution might be marginal for resolving relevant 
details of the penumbra, the initial condition in terms 
of a monolithic axisymmetric field is very artificial, and 
there is a strong influence from the (unavoidable) im- 
posed boundary conditions. Therefore, in this paper 
we investigate numerical resolution and boundary con- 
ditions. We also deepen our discussion on the flow field 
and magnetic structure in the observable layers of the 
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penumbra. At the bottom boundary most of the sunspot 
simulations impose a fixed magnetic field to prevent a 
rapid decay of the sunspot. Rempel ( 2011c| ) showed that 
this constraint can be relaxed in sufficiently deep do- 
mains where convective time scales reach several days, 
but not in shallow ~ 6 Mm deep domains typically used 
in simulations that address sunspot fine structure. Here 
we investigate how the magnetic top boundary condition 
as well as grid resolution influence details of penumbral 
fine structure, while we keep the initial state as well as 
the bottom boun dary unchanged . We deepen the anal- 
ysis presented in Rempel (2011b) to test the numerical 
robustness of their findings. In addition we expand the 
discussion to address three aspects that have been sub- 
ject to major controversy in recent years: 

(1) The convective nature of the penumbra: Vigorous 
convection (about half the strength of granulation) is 
the key process in MHD simulations leading to penum- 
bral fine structure. Indirect observational evidence for 
overturning motions is derived from 'twi sting' motions 
( Ichimoto et al.|2007b Bharti et al.| [2010) or correlation 
tracking that shows the convergence ot tr acers toward the 
edge o f filaments similar to g ranulation (Marquez et al. 



2006) 



tive ro 



Zakharov et al. 



( 2008 ) found evidence tor convec 
es in filaments oriented parallel to the solar limb 
and estimated that the observed velocities are sufficient 
to explain the penumbral brightness. Direct evidence 
for co nvective motions has been questioned by some au 



thors (Franz & Schlichenmaier||2009 Bellot Rubio et al. 



20101, while it was found by others in high resolution 



observations (Sanchez A lmeida et al. 2007; Joshi et al. 
2011; Scharme r^ al.|2011||s"charmer fc Henriques|201ip . 
| Joshi et al.| ( |2011[) found convective downnows mostly in 
the inner penumbra using the r ather deep forming CI 



5380 l ine. |Scharmer et al.| ( 2011 ); Scharmer & Henriques 
(20111 found consistent results using the CI 5380 and Fel 
6301 lines in terms of an anisotropic convection pattern 
with a correlation between intensity and vertical flow ve- 
locity not very different from quiet sun. Such correlations 
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were also previously indicated in the analysis of | Sanchez] 
Almeida et al] ( |2007 ). 

(Tj The presence of opposite polarity magnetic flux 
throughout the penumbra and its relation to downward di- 
rected mass flux: While inverse polarity flux is certainly 



prese nt in the outer penumbra (del Toro Iniesta et al. 
20011, it is debated if this is also the case in the inner 



penumbra. Sanchez Almeida (2005) found substantial 



amount of opposite polarity flux throughout the penum- 
bra, accounting to abo ut 58% of the uns igned magnetic 
flux of the penumbra, Langhans et al. (2005) did not 
find strong evidence in magnetogram data, possible due 
to resolution effects. Recently Franz (2011) found with 



Hinode data that about 40% of downfiows in the outer 
penumbra are in regions with opposite polarity flux and 
consider this a lower limit. 

(3) The magnetization of the Evershed flow: Spec- 
tropolarimetric observations point to a flow i n regions 
with a substantial magnetic field (see, e.g. Ilchimoto 



et al.|2008j[Borrero|2009 Borre ro k Solan ki 2008, pTU| 
see also discussions in Brummcll et al.l (|2008|); I'rhomai 



£ 



(2010). Some models suggest that penumbra 



nomas 
filaments 



are mostly field free 'ga ps' (Spruit & Scharmer 2006 - 
Scharmer fc Sprui"t| 2006 1 . A similar conclusion was also 



derived from an interpretati on of 'twisting filam ent mo- 
tions' as fluting instability (Spruit et al. 20101, leading 
to an upper bound for magnetic field in flow channels 
of about 300 G. MHD simulations however point toward 
strong 1 — 2 kG field in flow channels, in the outer penum- 
bra the field is even enhanced in flo w channels compared 
to the background ( Rempel|2011b I . 



2. NUMERICAL SETUP 

The numerical models were computed with an ex- 
panded version of the MURaM radiative MHD code 
( Vogler et al. 2005 1 . We introduced an artificial limi- 
tation of Alfven velocities to 60 km s _1 in order to pre- 
vent severe CFL time step constraints and implemented 
a new numerical diffusivi ty approach, both are described 

Rempel et al.| (|2009b[) . For the simulations presented 



here with grid sizes of up to 4.8 billion grid points we did 
additional performance enhancements including a major 
rewrite of IO and improved the scalability of the code 
to up to 24576 cores. Most of the simulations presented 
here were computed in the 1024 — 12288 core range. 

We present here a series of numerical sunspot mod- 
els which were all computed in a domain with a size 
of 49.152 x 49.152 x 6.144 Mm 3 . All models were ini- 
tialized with an axisymmetric self-similar magnetic field 
configuration described in Appendix A. For the mod- 
els presented here we used the parameters Bq — 6.4 kG, 
Rq = 8.2 Mm, and zq = 6.4185 Mm, leading to a sunspot 
with an initial flux of <f> = 1.2 ■ 10 22 Mx, a field strength 
of Bo = 6.4 kG at the bottom of our domain (z = 
Mm), dropping to 2.56 kG at the top [z = 6.144 Mm). 

While keeping this initial condition fixed we explore the 
dependence of the resulting sunspots on the top bound- 
ary condition and grid resolution. 

We use in the horizontal direction periodic boundaries, 
the bottom boundary condition is open for convective 
flows in regions with \B\ < 2.5 kG, but closed otherwise. 
The open boundary condition imposes a symmetric mass 
flux (all three components) in the ghost cells and extrap- 
olates the gas pressure such that its value at the bound- 



ary is fixed. In outflow regions the entropy is determined 
by the upstream values from within the computational 
domain, in inflow regions the entropy is set to a fixed 
value that leads to the correct solar energy flux under 
quiet Sun conditions. The closed boundary condition in 
regions with \B\ > 2.5 kG is implemented through an 
antisymmetric mass flux. The top boundary condition is 
closed (vertical mass flux antisymmetric) and stress free 
for horizontal motions. 

For the magnetic field we use a boundary condition 
that allows us to manipulate the inclination angle of the 
magnetic field. We focus on a single sunspot, which im- 
plies that the horizontal boundary imposes same polarity 
spots nearb y and one thus produc es only a very subdued 



penumbra (Rempel et al. 2009a). A relaxation of the 



horizontal periodicity is non-trivial and we thus decided 
to use a boundary that allows to control the inclination 
angle of the magnetic field to mimic different global field 
configurations, which is described in detail in Appendix 
B. A free parameter, a, gives for a = 1 a potential field 
extrapolation subject to horizontal periodicity, i.e. the 
field becomes vertical asymptotically. Values of a > 1 
lead to more horizontally inclined field, i.e. the hor- 
izontal field component is approximately enhanced by 
a factor of a compared to the potential field reference. 
The resulting field extrapolations are for a ^ 1 not force 
free outside the computational domain. They are not in- 
tended to give a realistic description for the global field 
topology above a sunspot, rather to explore the influence 
of different field geometries on the structure of a sunspot 
penumbra. 

All of the simulations presented here were computed 
with gray radiative transfer. To allow for a better com- 
parison with observations through forward modeling of 
spectral lines we computed 2 non-gray models that are 
derived from the gray simulations. We have a non-gray 
version of the simulation with 32 [16] km grid resolution 
that was advanced for 26 minutes and we have a non- 
gray model with 12 [8] km resolution that was restarted 
from the upper half of our 16 [12] km case and advanced 
for 15 minut es. The non-gray mo del with 32 [16] km was 
analyzed by |Bharti et al ( |2011 ) , more results will be 
presented in forthcoming publications. 

In the following analysis we will focus entirely on the 
gray simulations and discuss penumbral structure mostly 
through quantities extracted on constant t surfaces. In 
Figure [T] we display a snapshot from our highest reso- 
lution gray simulation (16 [12] km) that we consider in 
this paper. The simulation was computed with the a = 2 
boundary condition and evolved for 1 hour at the highest 
resolution. 



3. GLOBAL SUNSPOT PROPERTIES 

3.1. Influence from magnetic top boundary condition 

In Figure [2] we show 4 sunspot models computed with 
the boundary conditions a = 1,1.5,2,2.5. The case 
a = 1 (lower right quadrant) leads to a very subdued 
penumbra with a few isolated filaments, which is con- 
sistent with the penum bra structure that was found in 
Rempel et al. (2009a) in the direction where the peri- 
odicity imposes the same polarity sunspots. The other 
three cases show more developed penumbrae, the overall 
radial extent increases with the value of a. In Figure [3] 
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Figure 1. Top: bolometric intensity, bottom: magnetic field strength on a vertical cut through the center of the sunspot. Displayed is a 
snapshot from our highest resolution simulation (16 [12] km). An animation is provided with the online version. 



we show azimuthal averages of (a) intensity, (b) Evershed 
flow, (c) vertical and radial magnetic field, and (d) incli- 
nation angle. These quantities are in addition averaged 
for about 1 hour in time. Apart from the monotonic 
increase of the penumbra extent with a, it is also evi- 
dent that the most dramatic differences occur between 
the the potential field case (a — 1) and the other three 
cases with a > 1. For a = 1 the azimuthally averaged ra- 
dial flow velocity stays around 500 ms -1 , while all other 
cases have outflows with more than 3 kms -1 on average. 
Similarly the curves for the radial magnetic field as well 
as inclination angle do not differ as much between the 



a > 1 cases as compared to a — 1 and the rest. While 
the potential field case reaches only an average inclina- 
tion of 50°, all the other solutions with more extended 
penumbra and Evershed flow reach about 65°. 

This leads to the surprising conclusion that the poten- 
tial field case (which is the only physical, i.e. force free so- 
lution outside the computational domain) is a clear out- 
lier compared to the rest. This result has to be seen in the 
context of horizontal periodicity, which is not the proper 
horizontal boundary condition to describe the magnetic 
field above sunspots. A more reasonable (but compu- 
tationally more difficult to implement) boundary condi- 
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Figure 2. Influence of the magnetic boundary condition at the top of the domain on the radial extent of the penumbra. Left: bolometric 
intensity, right: magnetogram at r = 1. In each panel the 4 quadrants correspond to simulations performed with different magnetic top 
boundary conditions, (different values for a as described in Appendix B). The choice of a = 1 corresponds to a potential field extrapolation. 
Note that all simulations were performed in a 49 Mm wide domain, we show here only subsections. 
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Figure 3. Azimuthally averaged quantities at the r = 1 level for the different top boundary conditions shown in Figure [2] a) bolometric 
intensity, b) radial (Evershed) flow velocity, c) vertical (solid) and radial (dashed) field components, and d) field inclination angle. Increasing 
values for a lead to more extended penumbrae with faster Evershed flows, the most dramatic change is from a = 1 to a = 1.5. In terms of 
the inclination angle the potential field boundary falls short of the other solutions by about 10-15 degrees. 
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Figure 4. Influence of the numerical grid resolution on the properties of the penumbra. Left: bolometric intensity, right: magnctogram 
at r = 1. In each panel the 4 quadrants correspond to simulations performed with different resolution as indicated in the corners. Note 
that all simulations were performed in a 49 Mm wide domain, we show here only subsections. 
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Figure 5. Azimuthally averaged quantities at the r = 1 level for the different grid resolutions shown in Figure H] a) bolometric intensity, 
b) radial (Evcrshcd) flow velocity, c) vertical (solid) and radial (dashed) field components, and d) field inclination angle. 
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tion would be the free expansion of magnetic field into a 
half room (perhaps considering in addition spherical ge- 
ometry), which would lead naturally to a magnetic field 
with stronger horizontal field components. The bound- 
ary conditions with a > 1 emulate that behavior while 
maintaining horizontal periodicity. Note that simply in- 
creasing the horizontal extent of the domain for a = 1 
does not alleviate this problem, since the field remains 
asymptotically vertical. For the following investigation 
we will use a top boundary condition with a = 2. 

3.2. Influence from numerical grid resolution 

We modify the numerical grid resolution in the range 
from 96 [32] to 16 [12] km (horizontal [vertical] resolu- 
tion). The sunspot models computed with different grid 
resolutions are not fully independent from each other: 
The model with a resolution of 48 [24] km resolution was 
started from a snapshot of the 96 [32] km resolution sim- 
ulation evolved for 1 hour and ran for an additional 5.5 
hours. The model with a resolution of 32 [16] km resolu- 
tion was started from a snapshot of the 48 [24] km reso- 
lution run after 3.3 hours and ran for an additional 2.2 
hours. The highest resolution case with 16 [12] km resolu- 
tion was started from the last snapshot of the 32 [16] km 
case and evolved for an additional hour. Overall we com- 
pare sunspot models that have been evolved for about 6 
hours from the initial state. We did not evolve the higher 
resolution cases for the full length of time because of their 
computational expense (1 hour at 16 [12] km resolution, 
3072 x 3072 x 512 grid points, costs about 800,000 CPU 
hours on a CRAY XT-5). For the small length and asso- 
ciated short time scales in the photosphere about 1 hour 
at the highest resolution is sufficient to allow for the so- 
lution to adapt to the grid spacing. Structures in the 
deeper parts of the domain are already well resolved at 
lower resolution. 

Figure H shows the (gray) intensity and magnetogram 
at r — 1; Figure [5] shows the corresponding azimuthal 
averages of quantities at r = 1. Overall there is no sig- 
nificant influence from grid resolution on the radial ex- 
tent of the penumbra as well as the magnetic structure at 
t = 1, displayed in panels (c) and (d). Significant differ- 
ences occur with regard to the average intensity profile 
and Evershed flow shown in panels (a) and (b) . Only the 
higher resolution cases show the formation of a plateau 
like feature at around 0.7/0 in the inner penumbra. The 
average Evershed flow speed is increasing monotonically 
with resolution from a peak at around 2.3 kins -1 for 
96 [32] km resolution to about 4.6 kms -1 for 16 [12] km 
resolution. However, we see also here a clear sign of 
convergence: while the relative resolution changes from 
96 [32] to 48 [24] km and from 32 [16] to 16 [12] km are 
the same, the relative changes in the peak flow velocities 
are 1.55 and 1.1, respectively. The strong dependence of 
the Evershed flow velocity on grid resolution is a conse- 
quence from the driving mechanism that is conc entrated 
in a th in boundary layer just beneath r = 1 (Rempel 
2011b), which is difficult to resolve numerically (see also 
Section [7]). The umbra is with 0.3/0 rather bright for a 
sunspot with an umbral field stre ngth exceeding 3 kG. 
We return to this point in Section 8.3 



Figure [6] presents the penumbral fine structure at the 
r = 1 level for our highest resolution case (16 [12] k m res- 
olution). This figure is very similar to Figure 3 in Rem 



pel (2011b[), in which we displayed the fine structure at 



32 1 16 1 km resolution for the double sunspot simulation of 
Rempel et al. (2009a). The large degree of similarity un- 



derlines that the details of fine structure are not very sen- 
sitive to the numerical setup. Penumbral filaments show 
an enhancement of the radial magnetic field component 
(panel b), while the vertical field component is strongly 
reduced (panel c). In the outer penumbra we find a sub- 
stantial amount of opposite polarity flux, which will be 
characterized further in Section [4.31 The combination of 
a enhanced radial magnetic field component with a re- 
duced vertical field component results in the uncombed 
structure of the penumbra with a strong variation of the 
field inclination angle (panel d). Fast horizontal (Ev- 
ershed) outflows are found along the almost horizontal 
stretches of magnetic field (panel e). Overturning con- 
vection (panel f) is the underlying driving mechanism, 
which will be characterized further in Section [5] We pro- 
vide an animation of Figure [6 with the online version. 

Figure [7] presents the height dependence of the az- 
imuthally averaged Evershed flow and the vertical rms 
velocity. As already pointed out in Rempel (201 lb I, the 
Evershed flow peaks in the deep photosphere and falls off 
rapidly with height; we see the transition to an inverse 
Evershed flow at about r = 0.01. This figure is qualita- 



tively very similar to the observations reported by Bellot 
see Figure 8 therein) and by jBorrero 
igure Up we show the height depen 



Rubio et al. l (12006 
et all (12008 1.Tn - 



dence of the vertical rms velocity. We find essentially a 
drop by a factor of 2 between the r = 1 and r = 0.1 
levels, but no further drop toward higher layers. 

4.1. Robustness of penumbral fine structure 

To characterize fine structure further and quantify the 
resolution dependence of quantities we present in Figure 
M correlations between the Evershed flow velocity and 
(a) intensity (b) magnetic field strength as well as (c) 
radial and (d) vertical magnetic field components. All 
correlations are computed based on fluctuations around 
a smooth background that was obtained through a con- 
volution with a Gaussian having a FWHM of 1.5 Mm. 
At t = 1 the Evershed flow is found in the inner penum- 
bra in bright features, whereas the correlation vanishes 
toward the outer penu mbra (panel a). A similar behav- 
ior was also found by |Schlichenmaier et aT| (2005) and 



Ichimoto et al. 



( 2007ap . Trie correlation between Ever- 
shed flow and magnetic field strength is negative in the 
inner half, but positive in the outer half of the penum- 
bra. An Evershed flow in stronger magnetic fi eld regions 
in the outer pen umbra was also inferr ed by |Tritschler| 



4. PHOTOSPHERIC FINE STRUCTURE 



et al. (2007) and Ichimoto et al. (20081 from net circu- 
lar polarization (JNCP) variations with the viewing an- 
gle. The sign change in the Vr — \B\ correlation results 
from two effects: a positive correlation with the radial 
field component and a negative correlation with the ver- 
tical field component throughout the penumbra (Figure 
[8fc,d). While the negative correlation with the vertical 
neld dominates the inner penumbra, the positive correla- 
tion with the radial field component dominates the outer 
penumbra (see also Figure pi) . Except for the correlation 
between vr and B z that is with values about —0.7 quite 
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Figure 6. Sunspot fine structure at r = 1 for the highest resolution run (16 [12] km) with the a = 2 top boundary condition: a) bolometric 
intensity, b) radial field strength, c) vertical field strength, d) inclination angle, e) radial (Evershed) flow velocity, and f) vertical velocity. 
Black contours indicate in panel e) outflows with more than 10 kms -1 , and in panel f) supersonic downflows. The latter coincide mostly 
with strong inverse polarity patches in panel c). The dashed circles indicate R=8 and R=18 Mm, the dotted circle R=16 Mm. An animation 
is provided with the online version. 
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Figure 7. Height dependence of a) azimuthally averaged Evershed flow and b) vertical rms velocity. 
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Figure 9. Same correlations as shown in Figure[8] All quantities are evaluated at the r = 1 for different grid resolutions (color). 
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Figure 10. Bivariate probability density function for vr and |B| at r = 1 for the highest resolution case, 
penumbra (R = 10 to 13 Mm), panel b) shows the outer penumbra (_R = 13 to 16 Mm). 



Panel a) shows the inner 



significant all others reach in the deep photosphere only 
moderate levels of up to 0.5. All correlations decay very 
rapidly with height. 

To quantify the robustness of photospheric fine struc- 
ture we focus now on the correlations at r = 1 and study 
the resolution dependence in Figure [9] The lowest reso- 
lution case (96 [32] km) somewhat misses the vr — I cor- 
relation in the inner and the vr — \B\ correlation in the 
outer penumbra, all other grid resolutions produce com- 
parable levels of correlations. We can conclude that the 
magneto-convection process underlying pcnumbral fine 
structure is well captured starting from a resolution of 
about 48 [24] km. This resolution also produces an aver- 
age Evershed flow amplitude at this resolution that dif- 
fers less than 20% from our highest resolution case. As 
we will discuss in Section [7| the trend of increasing mag- 
netic field in flow channels (see vr — Br correlation) and 
the resolution dependence of the Evershed flow have a 
common origin. 

4.2. Magnetized flow channels 

The positive vr — Br correlation throughout the 
penumbra increases in amplitude with increasing reso- 
lution, which strongly points toward an active process 
supporting strong horizontal field in the flow channels as 
opposed to remnant magnetic field due to unavoidable 
numerical diffus ion effects as speculated by Nordl und fc:| 
Scharmer] ( |2010[ ). It was pointed out by |Rempel| Q201 lb| ) 
that the horizontal field originates mostly trom the hor- 
izontal shear in the sub-photospheric Evershed flow, i.e. 
the term B z 3 z vr in the induction equation (see Section 
[7]) . The robustness of this effect points toward a strongly 
magnetized Evershed flow in photospheric layers as it 
has been inferred from spectropolarimetric observ ations 



see. 



e.g. 



Ichimoto et al 



2010 Borrero'2DD0rT 



20081 |Borrero fc Solanki||2008l 

11s result is not necess arily in 
contradiction to th e 'ga ppy' penumbra model of |Spruit| 
& Scharmer ( 2006 1 and Scharmer & Spruit I d2006D since 



the strong horizontal magnetic held is confined to the 
thin boundary layer at r = 1 (Rempel 2011b). The 
value of \Br\ is in flow channels typically about a few 
100 G stronger than in the background. In Figure [TO] we 
present bivariate probability density function for vr and 
\B\ in the inner penumbra (R = 10 to 13 Mm) in panel 
a) and the outer penumbra (R = 13 to 16 Mm) in panel 



b). In the inner penumbra we find field strength of up 
to 4 kG in regions with low radial velocity, which corre- 
spond to the spines. Fast outflows have field strength of 
more than 2 kG, mostly due to Br. We see a similar be- 
havior in the outer penumbra with reduced overall field 
strength and less pronounced spines, but also here fast 
outflows are associated with field strengths around 2 kG. 
The extension toward low field strength is caused by a 
few granules present in the region 13 Mm < R < 16 Mm. 
In both, inner and outer penumbra, there is also a pop- 
ulation of inflows with substantial field strength. We see 
also a trend for increasing \B\ with increasing outflow 
velocities (vr > 0). 

4.3. Inverse polarity flux 

Inverse polarity magnetic flux embedded in the sunspot 
penumbra is an integral part of penumbral fine structure. 
Due to the substantial amount of overturning convection 
(see Sect. [5]) a certain amount of magnetic field lines has 
to turn back into the photosphere within the penumbra. 
Re turn flux along the outer penumbr a rim was first found 
by jWestendorp Plaza et al. (19971 and interpreted in 
the context of predictions tro m fl ux tube models such as 



Thomas fc Montesinos ( 1993 ) and Montesinos & Thomas 
~( 1997J). Fast (supersonic) downflows and r eturn flux in 



the interior of th e penumbra were found by del Toro Ini 
esta et al. ( 2001[ ). Substantial amounts of inverse polarity 
flux w erededuccd from observations by |Sanchez Almeida| 
( |2005[ ), which we will use as a reference for the analysis 
presented he re. This subject i s not w ithout controversy, 
for example |Langhans et al. ( 2005 ) did not find much 
evidence i n m agne tog rams at high spatial resolution. 

In Fi 
polarity 



Figs. [TTland 
ity flux founc 



12] we quantify the amount o f inverse 
Following Sanchez | 



in the penumbra. 
Almeida| ( |2005[ ) , we show in the left panels of both fig- 
ures the quantities R(\B Z \) (solid) and R{B Z ) (dashed) 
as function of radius, where (. . .) indicates the azimuthal 
average. The right panels show the ratio of the radially 
integrated signed and unsigned magnetic flux. Figure |11| 
shows how these quantities change with optical depth in 
the range from r = 1 to r = 0.001. Significant amounts 
of inverse polarity flux are only present in the deep photo- 
sphere, where we find integrated over the whole sunspot 
up to 11%, already at t = 0.1 only half of that remains. 
Figure [12] shows how these quantities vary with grid res- 
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Figure 11. Quantification of inverse polarity flux in the penumbra for di fferent height levels. Pa nel a) shows the quantities R.(B Z ) 
(dashed) and R(\B Z \) (solid) as function of radius (compare to Figure 13 in |Sanchez Almeida| ( |2005| >) . Panel b) shows the ratio of the 
radially integrated signed and unsigned fluxes. 
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Figure 12. Quantification of inverse polarity flux in the penumbra at r = 1 for different grid resolutions (color). Panel a) shows the 
quantities R(B Z ) (dashed) and R{\B Z \) (solid) as function of radius, panel b) the ratio of the radially integrated signed and unsigned fluxes. 
In panel b) the dotted (dashed) lines show for comparison results with different top boundary conditions a = 1.5(2.5) for the 48 [24] km 



olution as well as boundary conditions. Comparing the 
integrated values at R = 18 Mm in panel (b), we see an 
increase of the inverse polarity flux from abut 8% to 11% 
over the resolution range considered here (about half of 
the spread is due to our lowest resolution run). We also 
find that more extended penumbrae tend to have more 
opposite polarity flux than less extended ones (yellow 
dashed/dotted lines). Given the fact that the effective 
extent of the penumbra computed with a = 1.5 (2.5) is 
16 (19) Mm, the corresponding fractions of opposite po- 
larity fluxes are 6 (12)%. For a better co mparison with 
the results from Sanchez Almeida| (12005]), we present in 
Table [l] quantities that are comparable to those in their 
Table 1. Here the quantity $ z is defined analogous to 
theirs as 



$ z (i?i <R<R 2 ) = 2ir R(B z )dR 



(1) 



the signed flux is 50% of the unsigned flux. In contrast to 
Sanchez Almeida ( 2005 ) we do not find opposite polarity 
flux at photospheric levels in the umbra (although, some 
amount of opposite polarity field is present in umbral 
dots beneath the r = 1 surface in our simulations). 

The fact that the amount of return flux in simulations 
shows only little resolution dependence might come as a 
surprise, since most of this flux is present at small scales 
which vary somewhat with grid resolution. However, re- 
turn flux is a consequence of vigorous overturning convec- 
tion, which itself shows only little resolution dependence 
(see Section pi) . Overturning convection is a direct con- 
sequence of energy flux constraints that are the same in 
all cases considered regardless of resolution. 

4.4. Relation between opposite polarity regions and 
downflows 



Ri 



For computing the values in Table [T] we used R s = 18 
Mm for the radius of the sunspot. We present here val- 
ues computed for the highest resolution case. Our simu- 
lated sunspot is about a factor of 1.6 larger in terms of 
unsigned and a factor of 2 larger in terms of signed flux. 
While we find a large degre e of qualitative agreement 
with |Sanchez Almeida ( 2005 1 , the overall amount of in- 
verse polltfTtyTluxTalls short by a factor of about 3 (first 
line in table) . We find a substantial amount of opposite 
polarity flux in the region 12.7 Mm < R < 18 Mm, where 



Figure 13 1) presents the filling factors of regions with 
opposite polarity field B z < — 25 G (blue) and B z < 
-250 G (red). Both peak at about R = 16 Mm, which 
has been also indicated in Figure [6] through the dotted 
circle. The peak filling factors are 0.4 and 0.2 for these 
two thresholds. The filling factor of sup erso nic down- 
flows (green) reaches about 1%. In Figure fl3b) we eval- 
uate the downward directed mass flux in these regions 
relative to the total downward directed mass flux as func- 
tion of radius. Up to 60% (44%) of the returning mass 
flux is found in regions with B z < —25 G (B z < —250 
G), up to 5% in supersonic downflows. Integrated over 
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Table 1 

Quantification of inverse polarity flux 



Location 


Description 


Unsigned Signed 


® z (0<R<R s ) 


full spot 


$o = 1-18 x 10 22 Mx 0.89<J>o 


$ z (R a /2<R<R s ) 


penumbra 


$! = 0.44*o 0.76<J>i 


® z (Rs/V2<R<R s ) 


outer penumbra 


*2 = 0.19*o 0.49*2 


* z (0<iJ< Rs/2) 


umbra 


$3 = 0.54*o 1.00*3 
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Figure 13. a) Filling factors of regions with B z < — 25 G (blue), B z < —250 G (red) and supersonic downflows (green) as function of 
radius at r = 1. b) Fraction of downward directed mass flux found in these regions as function of radius. We show quantities computed 
for the resolutions of 16 [12] km (solid), 32 [16] km (dotted), and 48 [24] km (dashed). 

in the center of the penumbra values around — 20deg 
are typical. The magnetic field starts with an angle of 
35 deg at R = 8 and continues to point upwards to about 
R = 14.5 Mm. The minimum inclination the field reaches 
is —10 deg at R = 16 Mm. Except for the cores of penum- 
bral filaments (equivalent to upflow regions) in the inner 
most 2/3 of the penumbra, the inclination angles of flow 
and field are quite different: Flows in the penumbra are 
not just flows along the magnetic field, in particular not 
in downflow regions. To illustrate this point more clearly 
we present in Figure [l4|b) the average angle between flow 
and field given by arccos(|w • B|/(|iT||B|)). In the inner 
most 2/3 of the penumbra the angle between flow and 
field is 10 — 20 deg in upflow and 25 — 60 deg in downflow 
regions. These values are larger than those evident from 
panel a) since we measure here the average of the local 
misalignment, which is different from the misalignment 
of the average field and flow. 
The inclination o f flow and field agrees qualitatively 



the entire penumbra (from R — 8 to R = 18 Mm) about 
40% (27%) of the downward directed mass flux is found in 
regions with B z < —25 (B z < —250). Supersonic down- 
flows contribute only 2 — 3%. Solid, dotted and dashed 
lines indicate the respective quantities for the resolution 
levels of 16 [12], 32 [16], and 48 [24], respectively. 

At t = 1 the average velocity in the supersonic down- 
flow regions is —9.6 kms -1 , with fastest flows reach- 
ing — 15 kms . The average magnetic field strength 
is 2.7 kG. The average gas pressure is with 6.8 ■ 
10 4 dyne cm -2 about 30% lower than the typical photo- 
spheric pressure (at a constant height level the difference 
would be even larger), leading to an average intensity 
that is with 1.1/© clearly enhanced, in particular for a 
strong downflow. This has consequences for the I — v z 
correlation in the outer penumbra, which we will discuss 
further in Sections [5] and |6] 

The fact that most of the returning mass flux is not 
associated with opposite polarity magnetic flux indi- 
cates that it returns mostly by submerging magnetic field 
rather than by flowing along downward directed field- 
lines. To quantify this aspect further we present in Fig- 
ure [l4k) the inclination of flows and magnetic field rela- 
tive to the horizontal radial direction, i.e. the quantities 
7„ = arctan(w z / 'v r) and jb = &Tctan(B z / B r) . We av- 
erage these expressions in Evershed flow channels sepa- 
rated into upflow and downflow regions, i.e. we consider 
regions with vr > and v z > 250 or v z < —250 ms _1 . 
In upflow regions both field and flow show a very similar 
inclination dropping from about 30 deg at R = 8 Mm to 
10 - 15 deg at R = 15 Mm. For R > 15 Mm the inclina- 
tion for the magnetic field stays at around 10 deg, while 
the flow inclination increases to 40 deg. The situation is 
different in downflow regions. Values for the flow inclina- 
tion are —40 deg in the inner and outer penumbra, while 



with the findings of |Scharmer et al.| ( [2011 SOM). They 
found that in locally bright features (uptiows) inclina- 
tion angles of flow and field agree very well, while the 
magnetic field was close to horizontal or weakly upward 
p ointing in dow nflow regions (locally dark features) . 



Franz (2011 ) concluded from Hinode data that at least 
40% of penumbral downflows contain magnetic field with 
opposite polarity, similar to our finding. This agreement 
might be coincidental as a substantial fraction of up- and 
down- flowing mass flux still might be hidden in observa- 
tions. |Franz (2011) speculated that this is a lower limit 
and possibly all downflows have opposite polarity flux. 
We see only a moderate increase of this fraction with res- 
olution (for the resolutions of 48 [24], 32 [16], 16 [12] km 
we find 33%, 36%, 40%, respectively). In the inner 
penumbra mass is returning beneath the photosphere by 
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Figure 14. a) Inclination of magnetic field and flow with respect to radial direction in Evershed flow channels at r = 1, separated by 
upflow and dow nflow regions, b) Angle between v and B for the same regions. Different line styles indicate results from different resolutions 
as in Figure 1131 



submerging entire fieldlines (that still connect to the up- 
per boundary) along the lateral downflow lanes of fila- 
ments without requiri ng opposite pol arity flux to form 
(see also Figure 19 in Rempel 2011b). We see a mod- 
erate increase in this number with overall extent of the 
penumbra for a fixed resolution of 48 [24] km from 29% 
to 36% (comparing the a — 1.5 and a = 2.5 cases). If we 
extrapolate our results assuming that 100% of downflows 
have opposite polarity flux (for which we do not see any 
indicatio n here), we would obtain values similar to those 
found by Sanchez Almeida (2005). 

From our models we expect that more than half of the 
returning mass flux in the penumbra is found in regions 
that have the same polarity as the umbra of the sunspot. 
Examples for such downflow s were recently detected by 
Katsukawa & Jurcak (20101, although their connection 



to overturning motions is not evident from their obser- 
vations. 

4.5. Mass Flux associated with large scale flow 
component 

The mass flux in the penumbra can be formally sep- 
arated into a large scale mean flow with upflows in the 
inner and downflows in the outer penumbra and later- 
ally overturning motions. The former corresponds in the 
photosphere to the Evershed flow. Note that this separa- 
tion is only meaningful in terms of the overall mass flux 
balance. Individual flow elements have always a com- 
bination of radial outflow and lateral overturning mo- 
tions. To quanti fy both compon ents of the mass flux we 
define (see also Rempel 2011b) mf — (v z g > 0) and 
rn~ = (v z g < (J), where (. . .) denotes the azimuthal 
average. The unsigned vertical mass flux is then given 
through m to t = "m-t ~ rn 7 > while the locally unbalanced 
mass flux associated with radially separated regions of 
up and downflows (and corresponding horizontal flows 
in between) is given by m mcan = m[ + m~ . As a rela- 
tive measure for the mass flux contained in the azimuthal 
mean flow we consider 



^ R\' 



IdR 



r-R-2 



in: r\> 



IdR 



(2) 



Evaluating this expression on a constant height surface 
about 350 km beneath the r = 1 level in the plage re- 
gion (to make sure we stay below r = 1 in the inner 
penumbra) we obtain with R\ = 6 and R 2 — 18 Mm 



values of 15 — 16% for grid resolutions from 48 [24] km 
to 16 [12] km. Similarly we obtain values within this 
range for the two simulations at 48 [24] km resolution 
with different extent of the penumbra due to changes in 
the top boundary condition. Here we used a value of 
R 2 = 16 Mm and R 2 = 19 Mm for the a = 1.5 and 
a = 2.5 cases to consider the different extent of the re- 
gion occupied by the Evershed flow according to Figure|3j 
T his is also in agree ment with the value found previously 
by Rempel (2011b) for the double sunspot simulation of 
Rempel et al.|"( 2009a ) . 



In order to evaluate this quantity on a constant r sur- 
face we have to replace m z with the mass flux normal to 
r levels and include a geometric factor considering pro- 
jection effects due to inclined r surfaces. Both effects are 
considered by using the expression m T — m ■ Vf/|3 z f| 
instead of m z at the respective r level. Here f denotes 
a horizontally smoothed r, since we did not find a suffi- 
ciently balanced mass flux if we use r at maximum reso- 
lution. We compare here results that were obtained after 
convolving r in the horizontal direction with a Gaussian 
having a FWHM of 192 km. We find again compara- 
ble values for the resolution levels from 48 [24] km to 
16 [12] km. Values of e for the r levels of 1.0 (0.1) are 
14 - 16% (8 - 11%) (note that this is relative to the un- 
signed total mass flux at the respective level, which is at 
the r = 0.1 level about 25% of that at r = 1). The value 
at t = 1 is very close to what we found on a constant 
height surface about 300 km deeper, pointing toward ro- 
bustness of e with regard to the height (or r level) as 
well as grid resolution. 

5. OVERTURNING CONVECTION AND VISIBILITY OF 
CONVECTIVE SIGNATURES 

With the advent of magneto-convective models of 
sunspot penumbrae, the role of overturning convection 
in the penumbra has seen a controversial debate. 

On the theoretical side, overturning convection is the 
key process responsi ble for the energy transport and fila- 
mentation (see , e.g., Spruit fc Scharmer|[2 006 ; Scharmer 
|fc Spruit1|2 006 ; Hcinc mann et al.|12007[ |S charmcr et at 
20081 I Rempel et al.|2009a|b||KitiashviIi et al.|2009||Rem- 



pel||2011b| ) 

On the observational side, the evidence for overturning 
convection is controversial. While the existence of down- 
flo ws and upflows is evident (see, e.g., the direct evidence 
by |del Toro Iniesta et al. 2 001[ Franz fc Schlichcnmaier 
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Figure 15. Relation between azimuthally averaged bolometric intensity and vertical rms velocity (a) as well as rms intensity contrast 
(b) for different grid resolutions (color). In addition the dotted (dashed) lines show for comparison results with different top boundary 
conditions a = 1.5 (2.5) for the 48 [24] km case. In panel a) the black dashed line indicates a relation of the form / <x \Jv r z ma (r = 1). 
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Figure 16. a) Filling factors of regions with v z > 250 ms~ 1 (blue) and regions with v z < —250 ms —1 (red), b) Average upflow (blue) 
and downflow (red) velocity in those regions. All quantities are computed at the r = 1 level. Solid, dotted, and dashed lines show results 
from simulations with 16 [12], 32 [16], and 48 [24] km resolution, respectively. 

We find again only a very moderate dependence on the 
grid resolution as well as on penumbra extent (dashed 
and dotted yellow line) . Only the contrast of the lowest 
resolution case shows a significant deviations in regions 
with / < 0.7/0. If there is any systematic variation at 
all, it is a very moderate increase of the rms velocity with 
increasing grid resolution. Overall the predictions about 
the amplitude of overturning convection in the penumbra 
are very robust. Lower values for the amplitude of over- 
turning motions would be only possible in combination 
with a substantially higher contrast to satisfy energy flux 
constraints. 

While the overall amount of overturning convection 
as characterized by the vertical rms velocity is robust, 
there are differences in how much up- and downflows 
contribute. To characterize their contribution further 



|2009| ), their nature is debated: Do up- and downflows 
correspond to the radial endpoints of penumbral flow 
channels and are just the vertical component of the Ev- 
ershed flow, or do individual penumbral filaments show a 
flow pattern similar to that of granules with a substantial 
amount of mass flux turning over in lateral directions? 
With regard to the latter direct and indirect evidence 
was presented by Marqucz et al. (2006); Sanchez Almeida 
et al. (2007); Ichimoto ct al. (2007b); Kimmclc (2008); 
Zakha r ov etal (2008); Bharti ct al. (2010); Josh i et al.| 



2011 



2011 



); Scharmer et al.] ~i |201ip ; jScharmc r fe Henriques| 
Other investigations claim that these motions 



do not exis t ( Franz fc Schlichenmaier|2009 Bellot Rubio 



eraI1[20T0L 



Figure [T5| quantifies the vertical rms velocity and the 
relative rms intensity contrast for different grid resolu- 
tions. In panel (a) the dashed line indicates an approxi- 
mate relation of the form: 



Joe \/u™ s (r = 1) 



(3) 



which was already found by Rempel (2011b). This re- 
lation means that about 50% of the level of overturning 
convection found in the plage surrounding the sunspot is 
required to maintain an intensity of about 0.7/0 typical 
for the inner penumbra. This is possible since the rms 
contrast is with 25% about a factor of 1.5 larger than the 
contrast in the plage (17%). Note that the contrast val- 
ues are a few % higher in our gray simulations compared 
to a non-gray run. 



we show in Figure 16 i) the filling factors of regions 
with more than 250 ms _1 upflow (blue) and less than 
—250 ms _1 downflow velocity (red). Filling factors of 
up- and downflows are comparable in the plage region 
outside the sunspot with values of about 45%. The outer 
penumbra shows a clear preference for downflows with 
filling factor of 55%, while upflows occupy only 30% of 
the area. In the inner penumbra the upflow filling factor 
exceeds the downflow filling factor in the highest reso- 
lution cas e and is comparable for the other cases. In 
Figure 16 a) we show the respective average velocities in 
up- and downflow regions. In the plage region downflows 
dominate over upflows by a factor of about 1.5. This is 
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Figure 17. Average azimuthal extent of intensity and velocity features: regions with v z > 250 ms 1 (blue), regions with v z < —250 
ms -1 (red), and regions with intensity 0.05/0 above background (green). Panel a) shows the average, panel b) the area weighted average. 
Solid (dotted) lines correspond to the simulations with 16 [12] (32 [16]) km resolution. 
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Figure 18. Correlation between vertical velocity and bolometric intensity as function of a) height for a fixed resolution (16 [12] km) and 
b) resolution for a fixed height (r = 1). 



mostly an optical depth effect since we show quantities 
on the t — \ surface. Since up- and downflow velocity 
amplitude decline with height, the elevated r levels in up- 
flow regions sample lower velocities than the depressed r 
levels in downflow regions, while the difference on a con- 
stant height surface is typically less pronounced. This 
"apparent redshift" is not to be confused with the well 
known convective blueshift that originates in unresolved 
observations from the different intensity weighting of up- 
and downflow regions dependent on the the temperature 
sensitivity of the considered spectral line. In the inner 
penumbra upflow velocities dominate over downflows by 
a factor of up to 1.7 (despite the fact that also here the 
elevated r levels in upflow regions produce a bias toward 
downflows). For the currently explored resolution range 
we do not yet see a convergence of filling factors as well 
as mean up- and downflow velocities. It is likely that the 
trend toward upflow dominance in the inner penumbra 
will continue somewhat further with increasing resolu- 
tion. 

In Figure [17] we present the widths of intensity, up- 
and downflow features as function of radius. To this 
end we determine the average azimuthal extent of re- 
gions with intensity larger than 0.05/0 compared to the 
background that is defined through a smoothed intensity 
profile (Gaussian with a FWHM of 1.5 Mm). Up- and 
downflows are again defined through regions with more 
than 250 ms" 1 velocity amplitude. We consider here the 



following measures for the widths 
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where Si(r) are the detected widths at a given radial po- 
sition. For n — we obtain the regular average, while 
n = 2 can be considered as an "area weighted average" 
if we make the assum ptio n that the area of a feature 
scales with sf. Figure [l7| shows in panel a) (s(r))o and 
in panel b) (s(r-)) 2 . Regardless of the adopted measure 
we find that the widths of intensity and flow features 
are comparable in the plage region, while flow features 
are systematically smaller in the penumbra by a factor 
of up to 2. We do not see convergence of (s(r))o with 
grid resolution (we show here the 16 [12] and 32 [16] km 
cases by solid and dotted lines) in neither penumbra nor 
plage due to small scale features that dominate by num- 
bers in both regions. The area weighted average (s(r))2 
shows better convergence in particular for the width of 
intensity features. The most dramatic changes occur in 
flow features in the mid to outer penumbra. The latter 
indicates that even higher resolution than 16 km might 
be required to capture the width of flow features in the 
penumbra properly and that the current width of about 
100 (200) km has to be considered as an upper limit for 
the average (ar ea weighted average) width. 

In Figure [18| we present the v z — I correlation, which is 
very often used to characterize convective energy trans- 
port. The v z — I correlation reaches in the inner penum- 
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Figure 19. Bivariate probability density functions for intensity and vertical velocity computed for the simulation with 16 [12] km resolution. 
We denned four regions: The inner penumbra is from R = 10 to R = 14 Mm, the outer penumbra from R = 14 to R = 18 Mm, and 
plage is everything outside R = 20 Mm. In addition we computed the distribution function for plage regions with less than 500 G vertical 
field. The top row is computed with the vertical velocity at r = 1, the bottom row with the velocity at r = 0.1. We show the PDFs on a 
logarithmic scale in the range from 0.001 to 1 of the maximum value, contour lines show the levels of 0.1 and 0.01. Dashed horizontal and 
vertical lines show the average values for the considered sub regions. 



bra values similar to those found outside the sunspot 
(peak at around 0.6). Interestingly, the correlation van- 
ishes in the outer penumbra. We further find that a 
significant correlation is only present with the velocity 
at t = 1 inside the penumbra, whereas t — 0.1 leads also 
to a very moderate (0.3 — 0.4) outside the sunspot. We 
find here again only little resolution dependence in the 
range from 48 [24] to 16 [12] km. A v z — I correlation in 
the inner penumbra comparable to the plage surround- 
ing the sunspot and a very weak correlation in the outer 
penumbra are captured well starting from a resolution of 
48 [24] km. 

Note that the drop of the correlation in the outer 
penumbra is due to photospheric effects. If we compute 
instead the correlation between m z and (e; nt + p)/ Q + 
\/2v 2 (here e- m t and p denote the internal energy and 
gas pressure) on a constant height surface located about 
350 km beneath average t = 1 in the plage region, we find 
a constant value of 0.65 for R > 8 Mm. This indicates 
that rather moderate values of the / — v z correlation in 
the photosphere cannot be taken as an argument against 
convective energy transport, in particular not if strong 
magnetic field is present. 

6. PROBABILITY DENSITY FUNCTIONS OF 
PHOTOSPHERIC VELOCITY FIELD 

In Section [5] we discussed the correlation between in- 
tensity and vertical velocity. We found comparable val- 
ues (~ 0.6) in the inner penumbra and the strong plage 
surrounding the sunspot, while the correlation disap- 
pears in the outer penumbra. In Figure [19] we show for 
the highest resolution case (16 [12] km) bivariate proba- 
bility density functions for vertical velocity and inten- 
sity, computed (left to right) for the inner penumbra 



(R = 10 to R = 14 Mm), the outer penumbra (R = 14 
to R = 18 Mm), plage (R > 20 Mm) as well as plage 
regions with |i? 2 | < 500 G. We use here wider mask s for 
the inner and outer penumbra than in Section |4.2| since 
we want to include all of the fast downflows at the outer 
edge of the penumbra. The PDFs are computed for abso- 
lute intensity values (rather than fluctuations about the 
background) to allow for a direct comparison of intensity 
values between penumbra and plage. 

Comparing first the plage region and plage with \B Z \ < 
500 G, we see that the moderate correlation of 0.59 is 
mostly due to a population of very bright downflows 
that are associated with strong field. Masking out these 
regions increases the correlation to 0.73. The correla- 
tion for the inner penumbra has a value of 0.53 very 
close to that of plage, however, the shape of the PDF 
is substantially different. While the correlation in the 
plage regions results in about equal parts from bright 
upflows and dark downflows, it is in the inner penum- 
bra mostly due to bright upflows, some with brightness 
values reaching 1.4/q. Similarly to the plage region, the 
decorrelation in the outer penumbra is due to a strong 
population of bright downflows. In particular, most su- 
personic downflows have a brightness comparable to Iq. 
Most of these features have also strong opposite polar- 
ity magnetic field, indicating that the mechanism for the 
brightness enhancement is likely similar in penumbra and 
plage. 

The less pronounced tendency for downflows to be 
darker in the inner penumbra is in part due to the fact 
that the typical width of flow features is about half 
the width of intensity features according to Figure [TTJ 
While the upflows are mostly located centrally in bright 
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Figure 20. Velocity distribution functions in inner penumbra (from R = 10 to R = 14 Mm, blue), outer penumbra (from R = 14 to 
R = 18 Mm, green) and plage (R, > 20 Mm, red). Left to right we show the distribution functions for vertical, radial, and azimuthal flows; 
top to bottom the levels of r = 1, 0.1, 10 -2 , and 10 -3 . 

functions for penumbra and plage between their obser- 
vations and our simulations exist. They did find quite 
comparable results for the CI 5380 and Fel 6301 lines, 
while we see a rather strong height dependence for the 
probability density functions and resulting correlations 
between r = 1 and r = 0.1. The bivariate probability 
density functions do not show substantial differences be- 
tween the 48 [24] km and 16 [12] km cases in terms of the 
features discussed above. 

PDFs for v z , vr and v® are presented in Figure [20j 
At r = 1 the PDF for v z peaks in the plage region 
near v z = 1.5 kms -1 (the "typical" upflow velocity), 
while the PDFs for inner and outer penumbra peak near 
v z = kms . The PDF in the inner penumbra is 
skewed toward upflows, the PDF in the outer penumbra 
has a far extending tail of fast supersonic downflows. To- 



patches, downflows are located near the edge and are 
partially found in bright and dark features. 

Computing the PDFs based on the vertical velocity 
at t = 0.1 does not lead to any significant correlation 
in the penumbra, while the plage region shows values of 
0.38 and 0.55 excluding strong field regions. The core 
of the PDFs leads in most cases to a stronger correla- 
tion. This implies that smoothing the data increases the 
overall level of correlation found for both penumbra and 
plage. For example convolving the data by a Gaussian 
with FWHM of 150 km leads for r = 1 to correlations 
of 0.63, 0.19, 0.72, and 0.82. T hese values are compa- 
rable to those recently found by |Scharmer fc Henriques 
(2011|) for the interior penumbra and quiet sun. We see 



overall a good agreement with their results, although dif- 
ferences in the shape of the I — v z probability density 
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ward lower r values the shape of the PDF for all three re- 
gions becomes similar, while the PDF in inner and outer 
penumbra remains more narrow compared to the plage 
region. Note that the fast downflow wing in the outer 
penumbra extends toward higher layers than the fast up- 
flow wing in the inner penumbra, which is only present 
at t = 1 . A sim ilar behaviour was seen by |Scharmer fc| 



Henriques (2011) who found that downflows have a sim- 
ilar strength in CI 5380 and Fel 6301, while upflows are 
stronger in the deeper forming CI 5380 line. 

The PDF for vr is broader in the penumbra compared 
to plage at all r levels shown. For r = 1 and 0.1 it is 
skewed toward outflows, while the inflow part matches 
up with that of the plage region; toward r = 0.001 it is 
skewed toward inflows, while the outflow part matches up 
that of the plage region. In the inner penumbra the PDF 
shows a secondary peak at r = 1 at an outflow velocity 
of about 7—8 kms -1 , corresponding to the "typical" Ev- 
ershed flow velocity within flow channels (the azimuthal 
average is about a factor of two smaller) . 

At t = 1 the PDF for laterally overturning motions 
in outer penumbra and plage are similar, while the inner 
penumbra shows a more narrow PDF. In higher layers 
the PDFs for penumbra are more narrow than the PDF 
for plage. 

Observed distribution functions for the vertical veloc- 
ity in the range from —2 to 2 kms -1 were shown by 



Franz & Schlichenmaier ( 2009[ ) . They found in the quiet 
sun that upflows dominate over downflows for all velocity 
intervals. In the penumbra (they did not distinguish be- 
tween inner and outer penumbra) upflows dominate for 
\v z \ < 400 ms _1 , downflows for faster velocities. The 
overall width of the PDF in the penumbra was about 
half of the value found in the quiet sun. While our re- 
sults certainly show about a factor of 2 difference in the 
widths (see also the vertical rms velocity in Figure 15), 
other details might depend also on the utilized spectral 
line as well as resolution and require likely a comparison 
through forward modeling. 

7. RESOLUTION DEPENDENCE OF SUBSURFACE 
ENERGY EXCHANGE 

In Figure [21] we analyze subphotospheric processes re- 
sponsible forthe Evershed flow as well as the horizontal 
m agnetic field foun d in filament flow channels. Similarly 
to Rempel (2011b I we focus here on regions with v z > 



and vr > in the inner penumbra, where the strongest 
driving of the Evershed flow takes place. 

The left panels (a,d,g) present averages of magnetic 
field, velocity and temperature. Apart from the dramatic 
change of the Evershed flow speed we also find that the 
enhancement of Br at r = 1 is less pronounced at lower 
resolution, consistent with Figure [9J 

The middle panels (b,e,h) show energy conversion 
terms that play a central role in the m aintenance of the 
Evershed flow. As discussed in detail in Rempel (2011b I 
pressure buoyancy driving takes place in uptiow regions, 
while downflows are close to hydrostatic, which is oppo- 
site to the situation in non-magnetic convection. The 
pressure driving in upflows is diverted by the magnetic 
field into the horizontal direction where the acceleration 
takes place, while the pressure force itself does not have 
a strong component in the horizontal direction. To quan- 
tify the influence from grid resolution on this process we 



show the work by the vertical component of the pres- 
sure/buoyancy force (black), work by the vertical com- 
ponent of the Lorentz force (red), work by the radial 
component of the Lorentz force (blue) and radial com- 
ponent of the acceleration force (green). Dashed lines 
correspon d to simp l ified ex pressions that were discussed 
further in Rempel (2011b). The dotted line indicates 
the work by the radial component of the pressure force 
multiplied a factor of 10 for better visibility. Regardless 
of resolution we find the same balance, only the overall 
amplitudes of the various terms increase with increasing 
resolution. A quadratic increase with the peak Evershed 
flow velocity is expected, since the acceleration term is 
given in leading order by — (gvRV z d z VR) and v z itself does 
not change with resolution (tied to the fixed amount of 
energy transported). 

The panels on the right (c,f,i) show the contributions 
from the stretching term (red) advection term (black) 
and divergence term (blue) in the induction equation for 
the radial field component. Here it is most striking that 
the stretching term and in particular the contributions 
from B z O z vr (red, dashed) triples in amplitude ov er the 
resolution r ange explored. As discussed in detail in |Rem-| 
|pel| ( |2011b ), there is a strong connection between the Ev- 
ershed flow and the horizontal magnetic field in a thin 
boundary layer beneath the r = 1 surface. The subsur- 
face shear profile of the Evershed flow induces strong hor- 
izontal field, which is distributed along the r = 1 level by 
overturning convection and leads to a thin "outer shell" 
of the filament with strongly enhanced horizontal field. 
This boundary layer is also the place where strong mag- 
netic tension forces are in balance with vertical pres sur e 
and horizontal acceleration forces (see also Figure 22). 
This strong linkage between Br and vr is also evident 
from Figures [3] and [9j in which both the Evershed flow 
and horizontal magnetic field show the strongest reso- 
lution dependence with a trend toward faster flows and 
stronger field at higher resolution. 

8. DISCUSSION 

8.1. Role of top boundary condition 

We found a very strong influence of the magnetic top 
boundary condition on the extent of the penumbra. We 
do not obtain a penumbra at all using a potential field 
extrapolation since the horizontal periodicity implies an 
asymptotically v ertical field. This is in agreement with 
earlier results by Rempel et al.| (2009a), who only found 
extended penumbrae in between nearby opposite polar- 
ity spots. This rather puzzling result requires further 
investigation, we see three possible implications: 

1. This result is entirely due to the horizontal period- 
icity that is chosen for computational convenience. 

2. Other physical processes are still missing or a in- 
su fficiently resolved, like turbulent pumping as suggested 

by 

w 



Montesinos & Thomas 
rile processes along the 



1997) 



Brummell et al. 



(2008) 



ines of turbulent pum ping are 
certain ly present in our simulations (see, also Rempel 
2011b I and can be associated with the opposite polarity 
flux present in the outer penumbra, we did not find any 
indication of a strong resolution dependence here that 
would indicate that these processes could be underesti- 
mated. 
3. Hysteresis from our initial state could prevent 
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Figure 21. Horizontal averages between R=ll and R=14.5 Mm over regions with vr > and v z > as function of height following 
|Rempel] ( |2011b| |. The left panels (a,d,g) show Br (black, solid), B z (black, dashed), 1/2vji (blue, solid), v z (blue, dashed), and T (red). 
The center panels (b,e,h) show the energy conversion rates by different force components: pressure/buoyancy forces in vertical direction 
(black), Lorentz force in vertical direction (red), Lorentz force in horizontal direction (blue), and acceleration force in horizontal direction 
(green). Dashed lines refer to simplified expressions as explained in the text. The right panels (c,f,i) show different contributions in 
the radial induction equation: stretching (red), advection (black), divergence (blue), resistivity (green). The red dashed line shows the 
contribution from the vertical subsurface velocity shear {B z d z vji). The numerical grid resolution is increasing top to bottom: panels (a,b,c) 
96 [32] km, panels (d,e,f) 48 [24] km, and panels (g,h,i) 16 [12] km. 



the formation of a penumbra. To test this we con- 
ducted additional experiments with a sunspot model sim- 
ilar to our 48 [24] km resolution case, but in a larger 
73.728 x 73.728 x 9.2 16 Mm 3 domain (this model has 
been described also in Rempel ( 2011c[ )). After running 
that simulation for a total of 24 hours of solar time with 
our a = 2 boundary condition, we switched back to 
a = 1 (potential field). Within about 30 minutes the 
entire penumbra including the associated Evershed flow 
disappeared. After running for 2 hours with the poten- 
tial field boundary condition we switched back to the 
a = 2 case and the entire penumbra inclusive Evershed 
flow reappeared in about 30 minutes with properties es- 



sentially identically to the previous. This indicates that 
there is no hysteresis from within the convection zone 
in our numerical setup. This is not necessarily in con- 
tradiction to observed hysteresis effects in sunspots if 
the hysteresis would be caused by the global overlying 
magnetic field structure. Independent evidence against 
hysteresis from ou r initial state comes from flux emer- 
gence simulations (Cheung et al. 2010), which also fail 
in producing a penumbra if a potential field boundary is 
used in a periodic domain. Overall this indicates that 
the coronal magnetic field overlying sunspots has a po- 
tential feedba ck on penumbra l structure. This has been 
proposed by |Liu et al. (2005); Deng et al. (2005), who 
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observed changes in penumbral structure associated with 
solar flares. Our simulations certainly indicate that such 
feedback is possible and that in response parts of penum- 
brae could disappear or form on a time scale of about 
30 minutes. Additional support for a key role of the 
magnetic fie ld overlying the sunspot comes from[Shimizu 
et al. (2012), who discovered in a developing active re- 
gion a chromospheric precursor of penumbra formation. 
They observed the formation of a magnetic canopy in the 
chromosphere before the penumbra formation took place 
in the photosphere. The penumbra reached in the end an 
extent similar to the previously observed chromospheric 
signature. Their observation indicates that the magnetic 
field at chromospheric levels likely plays an active role in 
the penumbra formation process. 

8.2. Convergence with regard to numerical resolution 

Overall we find that penumbral fine structure is robust 
with regard to grid resolution over the explored range 
from 96 [32] to 16 [12] km. Even our lowest resolution 
case (which could be considered as an intentional ex- 
periment of under-resolving penumbral structure) gives 
still a good qualitative agreement. Quantitative agree- 
ment requires better than 48 [24] km resolution. This 
agreement is in terms of the structure and correlations 
of quantities in the penumbra in particular Evershed flow 
channels, opposite polarity magnetic flux and its relation 
to downflows, the level of overturning convection and 
the subsurface magneto convective processes responsible 
for the filamentation and driving of the Evershed flow. 
All resolution levels considered are not sufficient to re- 
solve turbulence within flow channels. However, we also 
find strong magnetization of Evershed flow channels as 
a robust result. Therefore it is doubtful that a substan- 
tial amount of turbulence would develop at photospheric 
levels. Studies of sunspot penumbrae at higher resolu- 
tion than considered here would be only meaningful if a 
nested or adaptive grid is used that allows for high lo- 
cal resolution, while the domain remains large enough to 
capture an entire sunspot. 

We do not use any explicit viscosity or magnetic re- 
sistivity in our simulations. The unavoidable hyper- 
diffusivity used is resolution and scale dependent. It 
scales at least linear with grid spacing near discontinu- 
ities or regions with monotonicity changes, but has higher 
order or is even switched off in well resolved regions. Our 
experiments on numerical resolution implicitly explore 
the dependence of the solutions on hyper-diffusivities. 
From power spectra we estimate that hyper-diffusion af- 
fects mostly features that are resolved by less than 6 grid 
points or have a scale of less than 100 km in our simula- 
tion with 16 [12] km resolution. This is consistent with 
our findings that details of up- and downflow structures 
that exist on these scales in the penumbra are not fully 
converged yet. 

8.3. Comment on diffusion 

We found in our simulations a rather bright umbra of 
0.3/0 for a sunspot with 3.3 kG central field strength 
(see Figure pi). Further investigation shows that the um- 
bral brightness is influenced to some degree by numerical 
diffusivities, in particular mass-diffusion. Using a modi- 
fied numerical approach with a lower mass-diffusion rate 



reduces the umbral brightness by 0.05 — 0.1 Iq. This dif- 
ference is mostly due to a different number density of um- 
bral dots as details of diffusivities influence the magnetic 
field threshold at which umbral dots become completely 
suppressed. We did not find a significant influence on 
other aspects of the solution. Numerical mass-diffusion 
mimics to some degree the effects of ambipolar diffusion, 
which indicates that some detail of umbral dots could be 
influenced by multi-fluid effects. 

8.4. Visibility of convective signatures in real 
observations 

The vertical rms velocity depends very weakly on grid 
resolution, showing a very moderate increase with in- 
creasing grid resolution. The average downflow velocity 
is larger than the upflow velocity in plage and the outer 
penumbra (by about a factor of 1.5), but in the inner 
penumbra upflows dominate over downflows by about a 
factor of 1.7. In the outer penumbra downflows have 
a substantially larger filling factor (55%) than upflows 
(30%), while the upflow filling factor dominates in the 
inner penumbra by a small margin. These differences be- 
tween up and downflow regions become more pronounced 
with higher grid resolution in the inner penumbra. The 
consequence of these asymmetries is that upflows (down- 
flows) tend to be easier to observe in the inner (outer) 
penumbra, i.e. with insufficient observational resolution 
one should expect to see a pattern of upflows in the inner 
and downflows in the outer penumbra. That raises the 
question of wh at a sufficient obser vational resolution is 
in this context. Bharti et al. (2011 ) used a non-gray ver- 
sion of the model in 32 [16] km grid resolution to address 
this question through forward modeling. They focused 
mostly on a few filaments in the inner penumbra and 
found that the upflow would be clearly visible at 0.14" 
observational resolution, while the visibility of the lateral 
downflows depends strongly on the spectral lines used. 
They found only 200 ms" 1 in the Fe 7090 line, while the 
CI 5380 line yielded up to 800 ms -1 downflow velocity. 
Based on the results presented here these numbers are 
likely still on the optimistic side, since the azimuthal ex- 
tent of flow structures is st ill d ecre asing with increasing 
grid resolution (see Figure 17) and Bharti et al. (2011) 
did not consider additional effects from stray light and 
potential line blends that complicate the situation in real 
observations. In that sense we do not see at this point a 
direct conflict with non-detectio n of overturning motions 
in penumbra filaments such as ( Franz & Schlichenmaier 



2009 Bellot Rubio et al. 2010 ~. Furthermore, recen t 



observations in the C I 5380 line by Joshi et al. (2011) 



Scharmer et al. (2011 ) found direct evidence tor such mo- 
tions, although their analysis depends strongly on stray 
light removal and complications due to the projection of 
the Evershed flow on the line of sight. 

Overall this indicates that observations are approach- 
ing the required resolution, which makes the outlook on 
the next generation solar telescopes such as NST, Gre- 
gor and ATST very promising. A recent example of new 
details that can be learned from comparing high resolu- 
tion observations with state of the art numerical models 



was presented by Steiner et al. (2010) using Sunrise ob- 
servations of the quiet sun. On the side of simulations 
the presence of overturning motions is robust as they are 
linked to very fundamental energy flux constraints. How- 
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ever, details still change, indicating that also here higher 
grid resolution will be required for an in depth compar- 
ison with observations through forward modeling. For 
the current simulations spectral line features that result 
from structures on scales of 100 km or less should be in- 
terpreted with caution, since details are likely influenced 
by numerical diffusion. 

9. CONCLUSIONS 



In deepening our previou s studies in Rempel et al. 
( 2009b|a ) ; Rempel| ( 2011a|b ) we investigated how bound- 
ary conditions and grid resolution affect the structure of 



ary conditions and grid resolution attect the structure of 
penumbrae in simulated sunspots and worked out obser- 
vational consequences that can help to further constrain 
the magneto-convective structure of the penumbra. The 
main findings are the following: 

• The magnetic top boundary condition appears to 
determine the radial extent of the penumbra. We 
do not get an extended penumbra if we impose 
a potential field subject to horizontal periodicity. 
Thus, simulations in periodic domains cannot pre- 
dict the extent of penumbrae from first principles 
at this point. 

• Magneto convection forming penumbral fine struc- 
ture and driving the Evershed flow is robust over 
the explored resolution range. While all ex- 
plored grid resolutions reproduce basic aspects of 
a penumbra, a resolution of at least 48 [24] km 
(horizontal [vertical]) is required for a quantitative 
agreement. 

• The overall amount of convective motions in the 
photosphere as characterized by the vertical rms 
velocity shows little dependence on resolution and 
extent of the penumbra. It is linked to the intensity 
through an approximate relation of the form / oc 
^"(t = 1). 

• We find asymmetries between upflow/downflow fill- 
ing factors and velocities that favor the detection 
of upflows in the inner and downflows in the outer 
penumbra. These asymmetries become more pro- 
nounced at higher resolution. Forward modeling 
likely overestimates the visibility of a convective 
flow pattern in penumbrae at this point. 

• Evershed flow channels are strongly magnetized 
at photospheric levels, values around 1.5-2.5 kG 
are typical. The tendency for enhanced horizontal 
magnetic field within flow channels by a few 100 G 
compared to the background is a robust feature. 

• Opposite polarity flux is present in the penumbra 
at a level of about 10% of the total flux of the 
sunspot, with no resolution dependence for an oth- 
erwise fixed setup. The overall amount of opposite 
polarity flux increases with the extent of a penum- 
bra from 6 — 12%. 

• Integrated over the entire penumbra about 40% of 
the downward directed mass flux is associated with 
magnetic flux having a polarity opposite to that of 
the umbra. We see only a moderate increase of this 



number with resolution and extent of the penum- 
bra. Thus, there should be a significant fraction 
of downflows in the penumbra with the same field 
polarity as the umbra, i.e. a substantial amount of 
mass flux returns by subduction of upward pointing 
magnetic field rather than by flowing along down- 
ward pointing field lines. 

• While strong supersonic downflows are present in 
the outer penumbra (we found up to 15 kms -1 flow 
amplitude) , their overall contribution to the down- 
ward directed mass flux is with 2 — 3% negligible. 

• In the outer penumbra we find a substantial popu- 
lation of bright downflows that lead to a vanishing 
I — v z correlation, while the correlation in the inner 
penumbra and surrounding plage region are compa- 
rable (reaching moderate values around 50 — 60%). 
Nevertheless, energy is transported in all these re- 
gions by convection. Thus, low or moderate I — v z 
correlations cannot be taken as an indication for 
the absence of convective energy transport, in par- 
ticular when magnetic field is present in the pho- 
tosphere. 

• Most of the vertical mass flux in the penumbra is 
balanced within each radial shell. The unbalanced 
component of the mass flux (upflows in the inner 
and downflows in the outer penumbra) is about 
15% of the total overturning mass flux integrated 
over the whole penumbra. 

We conclude this paper by summarizing the magnetic 
field and flow structure of penumbral fil ame nts found in 
numerical simulations to date in Figure [22j This sketch 
is mos tly based on the conclusions from this paper as 
well as Rempel ( 2011b| ) (see Figure 17 therein). 

The side view shows the central upflow regions of fila- 
ments where the average field and average flow are well 
aligned. The pressure driving in the upflows is almost 
orthogonal to the mostly horizontal acceleration of the 
Evershed flow, requiring the Lorentz force to facilitate 
the energy exchange (the net work by the Lorentz force 
is essentially zero). In the lateral downflow regions the 
magnetic field continues to point upward leading to a 
substantial misalignment between flow and field and sub- 
mergence of field lines. The side view describes mostly 
the situation in the inner penumbra, where we do not 
find a substantial amount of inverse polarity magnetic 
flux. 

The top view shows the central upflow region with the 
adjacent lateral downflows. While mass is moving out- 
ward (Evershed flow) it moves to the edge of filaments 
where it returns beneath the solar surface mostly through 
submergence of field lines. Only near the outer end of the 
filament a significant fraction of the mass flux returns 
through flows along downward pointing field lines. 

The cross section shows a structure with reduced mag- 
netic field strength in the subphotospheric layers that is 
filled with an overturning convection pattern maintain- 
ing the penumbral brightness. While this structure has a 
lot of similarity with the i d ealized field free gap mode l of 
Spruit fcScharmer| ( [2006| ; |Scharmer fc Spruit| ( |2006[ ), it 
is not field tree. The filament has a core with a su bstan- 
tial vertical field component (see also Figure 17 of Rem- 
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Side view 




Bx, Bz 



a i, 



Penumbral filament in MHD simulation 





Top view 



Cross Section 



y the direction along the filament away from the spot center. In the central upfiow regions flow and ficfd arc weff aligned, while the flow 



Figure " 22. A sk etch summarizing the basic field and flow structure of a penumbral filament as present in the numerical simulation. We 
present a schematic side, top and cross-section view, x and z denote the horizontal and verti cal direction perpendicular to the filament 



submerges mostly horizontal held lines in lateral downtlow regions. Overall filaments have a reduced field strength, but they contain a core 
with a non-vanishing vertical field component. Some ol the associated Hux continues upward, some of the Hux returns downward within the 
filament cavity. Depending on the position of the r = 1 fevel the latter might become visible as inverse polarity tlux. The strong subsurface 
shear of the Evcrshcd flow induces a strong horizontal field component that is concentrated along the r = 1 surface. This leads to strongly 



magnetized Evcrshcd How channeIs"Tr7TrTevisiblc layel^^hlle'TneTreTorstrength is significantly rcducecTTrTTrTe'subsurlacc layers. 



pel (2011b)) that is spreading out laterally near the top 
ol the tilament. A fraction of it passes through the curved 
r = 1 surface (indicated by a dashed line) and contin- 
ues upward. Overturning convection can bend some of 
these field lines back downward, leading to opposite po- 
larity flux along the side of filaments. If the r = 1 level 
is located deep enough, this can lead to the appearance 
of opposite polarity flux along the edge of filaments in 
magnetograms. The Evcrshcd flow (shown in green) in- 
creases rapidly toward the r = 1 level in a region with a 
substantial vertical field component, leading to induction 
of a strong horizontal field component along the axis of 
the filament. Overturning motions distribute the mag- 
netic field along the r = 1 level in a thin boundary layer 
shown in red, resulting in strongly magnetized Evershed 
flow channels at photospheric levels that might appear 
like horizontal flux tubes to the observer. Note that we 
simplified the picture in the sketch with regard to the 
distribution of the horizontal field component. Horizon- 
tal field is present throughout most of the filament, it 
is however about a factor of 2 stronger in the boundary 
layer indicated in the sketch and mostly there dynami- 
cally relevant (the acceleration of the Evershed flow takes 
place in this boundary layer, see indicated force balance). 
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10. APPENDIX A 



For initialization of our simulations we use an axisymmetric self-sim ilar magnetic field configuratio n . This ap pr 
has been used previo usly for magnetohydrostati c sunspot models by Schluter & Temesvary (19581; Deinzer| |Tl 



and more recently by Schussler & Rempel (2005). 



roach 
1965), 



B z (R,z)=B f(0g(z), 
B R (B,,z) = -Bo^f(()g'(z) 



(5) 
(6) 



Here R, z correspond to the radial and vertical direction in cylindrical coordinates, g(z) is the vertical magnetic field 
profile along the symmetry axis and £ = R^g(z). The function /(C) describes the height independent radial profile of 
the vertical field component and can be freely chosen. For initializing our simulations we used the following functions 
to describe the radial and vertical field profile: 



/(C) = cxp[-(C/ J Ro) 4 ] 
g(z)=exp[-(z/z Q ) 2 }. 



(7) 
(8) 



leading to a total of three parameters that determine the field structure: Bq, i?o, and z$. The total flux is given by 
$o = 2.7833-Bo-Ro- Note that the most important parameters for our problem are $o and Bo, the detailed choices 
of /(C) and g(z) are secondary. We have chosen g(z) in this form to ensure g'(z = 0) = for consistency with our 
vertical field bottom boundary condition and /(C) to initially concentrate most field into the sunspot. The typically 
used Gaussian profile would lead to a smaller sunspot and a stronger plage region. 

We insert the above magnetic field structure into a thermally relaxed HD simulation. We do not correct the thermal 
structure to account for magnetic forces, instead we allow this adjustment to happen as part of the the dynamical 
evolution. We set however in regions with B > 2 kG the entropy to the value found at the bottom of the domain, 
while keeping the pressure unchanged and initially quench velocities in magnetized regions. The latter two have a 
rather minor effect on the overall evolution. 
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11. APPENDIX B 

For the simulations presented in this publication we used a modified top boundary condition, which allows us to 
arbitrarily to change the inclination angle of the magnetic field, while maintaining horizontal periodicity. To explain 
this boundary condition we start fro m a potential fi eld extrapolation in a horizontally periodic domain, which can be 
derived from a potential of the form ( Cheung||2006 ): 

$(z, y,z)= f ${k x , k y )e-^ z e lk * x e lk y y dk x dk y , (9) 



with \k\ = \/k^ + k 2 The magnetic field components are given by B — — V<&, together with the relation Bo(k x , k y ) = 
|fc|$(fca;, ky). Here B$ is the Fourier transform of B z (x,y,0) at the boundary z — 0. 

B * = -J Ba(k„ k v ) l -^re-^ z e^ x e^ y dk x dk y 
B y =- J B (k x , k y )^e-^ z e ik ^e ik ^dk x dk y 

B z = + J B Q (k x , k v )e-\ k \ z e ik * x e ik *vdk x dk y (10) 

From this potential field extrapolation we derive a generalized non-potential field boundary condition by introducing 
an additional parameter a controlling the field inclination angle: 

B« = -aJ B (k x , k y )^e- a ^ z e^ x e tk y y dk x dk y 
B v=- a f £o(*b, k y )^e- a W*e ik ' x e ik * v dk x dk v 

B* = + f B Q (k x ,k y )e~ aWz e tk ^e ik " v dk x dk y (11) 

While using a value of a = 1 reproduces the potential field, a value of < a < 1 leads to a field more vertical than 
potential (vertical field for a — 0), a value of a > 1 to a field extrapolation more horizontal than potential. 
The electrical current distribution associated with this field extrapolation is given by: 

j" = +(l -a 2 ) I B {k Xl k y )ik y e- aWz e lk * x e ik y y dk x dk y 

3 y = -(1 - a 2 ) J B {k x , k y )ik x e- a W z e lk * x e tk y y dk x dk y 

tf = (12) 

While the vertical current remains zero in the region z > 0, closed horizontal current loops induce magnetic field 
which changes the effective field inclination in z < 0. The magnetic field in the region z > is not force free, but the 
magnetic field in the computational domain z < remains close to force free due to the low ft conditions we consider. 



